{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Design a 15th-order bandpass filter in SOS format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import signal\n",
    "import numpy as np\n",
    "#sos = signal.ellip(15, 0.5, 60, (0.2, 0.4), btype='bandpass', output='sos')\n",
    "sos = signal.cheby1(15, 0.5, (0.2, 0.4), btype='bandpass', output='sos')\n",
    "#sos = signal.cheby2(15, 60, (0.2, 0.4), btype='bandpass', output='sos')\n",
    "# sos = signal.bessel(15, (0.2, 0.4), btype='bandpass', output='sos')\n",
    "# sos = signal.butter(15, (0.2, 0.4), btype='bandpass', output='sos')\n",
    "sos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the frequency response at 1500 points from DC to Nyquist."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w, h = signal.sosfreqz(sos, worN=1500)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the response."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.subplot(2, 1, 1)\n",
    "db = 20*np.log10(np.maximum(np.abs(h), 1e-5))\n",
    "plt.plot(w/np.pi, db)\n",
    "plt.ylim(-75, 5)\n",
    "plt.grid(True)\n",
    "plt.yticks([0, -20, -40, -60])\n",
    "plt.ylabel('Gain [dB]')\n",
    "plt.title('Frequency Response')\n",
    "plt.subplot(2, 1, 2)\n",
    "plt.plot(w/np.pi, np.angle(h))\n",
    "plt.grid(True)\n",
    "plt.yticks([-np.pi, -0.5*np.pi, 0, 0.5*np.pi, np.pi],\n",
    "           [r'$-\\pi$', r'$-\\pi/2$', '0', r'$\\pi/2$', r'$\\pi$'])\n",
    "plt.ylabel('Phase [rad]')\n",
    "plt.xlabel('Normalized frequency (1.0 = Nyquist)')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the same filter is implemented as a single transfer function, numerical error corrupts the frequency response:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#b, a = signal.ellip(15, 0.5, 60, (0.2, 0.4), btype='bandpass', output='ba')\n",
    "b, a = signal.cheby1(15, 0.5, (0.2, 0.4), btype='bandpass', output='ba')\n",
    "#b, a = signal.cheby2(15, 60, (0.2, 0.4), btype='bandpass', output='ba')\n",
    "# b, a = signal.bessel(15, (0.2, 0.4), btype='bandpass', output='ba')\n",
    "#b, a = signal.butter(15, (0.2, 0.4), btype='bandpass', output='ba')\n",
    "w, h = signal.freqz(b, a, worN=1500)\n",
    "plt.subplot(2, 1, 1)\n",
    "db = 20*np.log10(np.maximum(np.abs(h), 1e-5))\n",
    "plt.plot(w/np.pi, db)\n",
    "plt.ylim(-75, 5)\n",
    "plt.grid(True)\n",
    "plt.yticks([0, -20, -40, -60])\n",
    "plt.ylabel('Gain [dB]')\n",
    "plt.title('Frequency Response')\n",
    "plt.subplot(2, 1, 2)\n",
    "plt.plot(w/np.pi, np.angle(h))\n",
    "plt.grid(True)\n",
    "plt.yticks([-np.pi, -0.5*np.pi, 0, 0.5*np.pi, np.pi],\n",
    "           [r'$-\\pi$', r'$-\\pi/2$', '0', r'$\\pi/2$', r'$\\pi$'])\n",
    "plt.ylabel('Phase [rad]')\n",
    "plt.xlabel('Normalized frequency (1.0 = Nyquist)')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
